%ErrorCheck.m
function ErrorCheck(orig,Tanalog_Original,rec,Tanalog_rec,fmax,OrigTimeAxis,SNR,Type,H2)
% orig = Original_Signal;
% rec  = real(CurRunDB.RecSignal);

%For Plotting use PrintFlag=1;
PrintFlag=1;
% NOTE:
% for running offline use the next command (after running start.m the first time):
% ErrorCheck(Original_Signal,Tanalog_Original,real(CurRunDB.RecSignal),1/f_rec,fmax,OrigTimeAxis,SNR,Type);

close all;

FsOrig=1/Tanalog_Original;
FsRec=1/Tanalog_rec;

[orig_vec,orig_freq] = FreqDomain( FsOrig , orig, fmax);
[rec_vec,rec_freq]   = FreqDomain( FsRec, rec ,fmax);

%Calibration Method
CurrentCosFreqStruct = load('CurrentCosFreq.mat');
CurrentCosFreq = CurrentCosFreqStruct.CurrentCosFreq;
H_2=0;
H2 = CalibrationFunction(H_2,CurrentCosFreq,orig_freq,rec_freq);

%% Frequency Domain:

rec_freq = interpft(rec_freq,length(orig_freq));
rec_vec = interpft(rec_vec,length(rec_freq));

% Smooth signals
smooth_factor = 3;
orig_freq=smooth(abs(orig_freq),smooth_factor);
rec_freq=smooth(abs(rec_freq),smooth_factor);

% Amplitude zero order calibrating: %NOTE: ask Debby
max_orig=max(orig_freq);
max_rec=max(rec_freq);
rec_freq=rec_freq*max_orig/max_rec;

% Plotting signals and diff in freq domain
if (PrintFlag)
    figure;
    subplot(2,1,1);
    % % dB plotting:
    plot(rec_vec,db(abs(rec_freq)),orig_vec,db(abs(orig_freq)));
    % % linear plotting:
%     plot(rec_vec,(abs(rec_freq)),orig_vec,(abs(orig_freq)));
    title('Recovery signal and Original signal - Freq domain');
    legend('Recovery Signal','Original Signal');
    xlabel('Frequency (Hz)');
    ylabel('Amplitude [dB]');
    xlim([-fmax,fmax]);
    % Ploting diff between Rec and Orig
    diff = abs(rec_freq)-abs(orig_freq);
    % error = norm(diff);
    subplot(2,1,2);
    plot(orig_vec,diff);
    title('Diff = |Recovery signal| - |Original signal| - Freq domain');
    legend('Diff = |Recovery signal| - |Original signal| ');
    xlabel('Frequency (Hz)');
    ylabel('Amplitude');
    xlim([-fmax,fmax]);
end
%% Time Domain: (New way of determining time domain)

% plotting the 2 signals before shift:
% figure;
% subplot(2,2,1:2)
% plot(orig);
% axis tight;
% subplot(2,2,3:4)
% plot(rec);
% axis tight;

[P1,Q1] = rat(FsOrig/FsRec);% Rational fraction approximation
FsOrig/FsRec
rec = resample(rec,P1,Q1);        % Change sampling rate by rational factor

% Finding a Signal in a Measurement
% We can now cross-correlate signal S to templates T1 and T2 with the xcorr function to determine if there is a match.

[Corr,lag] = xcorr(rec,orig);

if (PrintFlag)
    figure;
    subplot(2,2,1:2);
    plot(lag/FsOrig,Corr,'k'); ylabel('Amplitude'); grid on
    title('Cross-correlation between Rec and Signal - before cyclic shift')
    xlabel('Time(secs)');
    [~,I] = max(abs(Corr));
    TimeDiff = lag(I)/FsOrig;
    Tindex = lag(I);
    cprintf('blue','Signal type: %s\n',Type);
    fprintf('lengths before paddind: orig: %d   rec: %d\n',length(orig),length(rec));
    % shifting the signal to the correct place
    % Option 1: cyclic shift (not so good)
    % rec = circshift(rec.',-Tindex).';
    % Option 2 - (Better) zerro padding
    if (Tindex < 0)
        fprintf('Tindex = %d < 0\n',Tindex);
%             rec = [zeros(1,abs(Tindex)),rec];
%             orig(1:abs(Tindex))=0;
        if ((length(rec)-Tindex) > length(orig)) % for notoverflowing the array
            FinalIndex=length(orig);
        else
            FinalIndex=length(rec)-Tindex;
        end
        orig = orig(-Tindex+1:FinalIndex);
        
    else
        %     fprintf('Tindex = %d > 0\n',Tindex);
            rec = [rec,zeros(1,Tindex)];
            orig= [orig(1:length(rec)-Tindex), zeros(abs(Tindex))];
        
    end
    fprintf('lengths after paddind: orig: %d   rec: %d\n',length(orig),length(rec));
    % plotting Corr after shift:
    [Corr,lag] = xcorr(rec,orig);
    subplot(2,2,3:4);
    plot(lag/FsOrig,Corr,'k'); ylabel('Amplitude'); grid on
    title('Cross-correlation between Rec and Signal - after shift')
    xlabel('Time(secs)');
end

%% Calculating correlation
% This notion of similarity is similar to rote correlation of the 
% two sequences but is invariant to time delay.

norm_max_xcorr_mag = @(x,y)(max(abs(xcorr(x,y)))/(norm(x,2)*norm(y,2)));

%%
% Example:
% x=[1,2,3,7,9];
% y=[0,1,2,3,5];
% norm_max_xcorr_mag(x,y);
if (SNR)
    fprintf('SNR = %d : \n',SNR);
end

% In Freq Domain:
x = orig_freq; y = rec_freq;
ErrFreq=norm_max_xcorr_mag(x,y);
cprintf('blue','correlation in Freq domain: %.3g    (smooth_factor = %d)\n',ErrFreq,smooth_factor);

% In time domain- see below
    % old way:
    % x = downsample(orig,floor(length(orig)/length(rec))); y = rec;
    % ErrTime=norm_max_xcorr_mag(x,y);
    % fprintf('correlation in Time domain: %.3g\n',ErrTime);

% Amplitude zero order calibrating: %NOTE: ask Debby
max_orig=max(orig);
max_rec=max(rec);
rec=rec*max_orig/max_rec;

if (PrintFlag)
    % plotting the 2 signals:
    figure;
    subplot(3,1,1);
    plot(orig,'g');
    title('Original signal - Time domain');
    legend('Original Signal');
    xlabel('time (sec)');
    ylabel('Amplitude');
    axis 'tight';
    % hold all; hold on;
    
    subplot(3,1,2);
    plot(rec,'b');
    title('Recovery signal - Time domain - after shift');
    legend('Recovery signal - after shift');
    xlabel('time (sec)');
    ylabel('Amplitude');
    axis 'tight';
    
    
    
    subplot(3,1,3);
    plot(rec,'b');
    hold on;
    plot(orig,'g:');
    title('Recovery signal and Original signal - Time domain - After shifting');
    legend('Recovery Signal','Original Signal');
    xlabel('time (sec)');
    ylabel('Amplitude');
    axis 'tight';
end


% Now calculatin Correlation In time domain:
x = orig; y = rec;
ErrTime=norm_max_xcorr_mag(x,y);
cprintf('blue','correlation in Time domain: %.3g\n',ErrTime);


%% Play sound of signals - %% ETGAR - takes time and quite loud. wait patiently 
% wavplay(orig);
% pause(2); % wait for 2 seconds quietly before playing the rec_signal;
% wavplay(rec);


%% now we examine the coorelation between Ofiginal and a random signal:
% fprintf('\nComparing with a RANDOM SIGNAL:\nRatio is the Ratio between the correlation aforementioned and the correlation with the Random signal\n\n');
% max_amp=max(orig);
% rand_seq = randsrc(1,length(orig),linspace(-1,1,1000));
% random_rec = (max_amp/2)*rand_seq;
% 
% % Time domain:
% x = orig; y = random_rec;
% ErrTimeRand=norm_max_xcorr_mag(x,y);
% fprintf('correlation between original and random signal - Time domain: %.3g\n',ErrTimeRand);
% if (ceil(ErrTime/ErrTimeRand) >= 10)
%     TimeOrderOfMagnitude = ceil(log10(ErrTime/ErrTimeRand))-1;
%     cprintf('blue','In Time Domain: Ratio of ~%d (more than %d order of magnitude)\n\n',ceil(ErrTime/ErrTimeRand),TimeOrderOfMagnitude);
% else
%     cprintf('red','WARNING: In Time Domain: Ratio of ~%d --> less than 1 order of magnitude!\n\n',ceil(ErrTime/ErrTimeRand));
% end 
% 
% % Freq Domain:
% [random_rec_vec,random_rec_freq]   = FreqDomain( 1/Tanalog_Original, random_rec ,fmax);
% x = orig; y = random_rec_freq;
% ErrFreqRand=norm_max_xcorr_mag(x,y);
% fprintf('correlation between original and random signal - Freq domain: %.3g\n',ErrFreqRand);
% if (ceil(ErrFreq/ErrFreqRand) >= 10)
%     FreqOrderOfMagnitude = ceil(log10(ErrFreq/ErrFreqRand))-1;
%     cprintf('blue','In Freq Domain: Ratio of ~%d (more than %d order of magnitude)\n\n',ceil(ErrFreq/ErrFreqRand),FreqOrderOfMagnitude);
% else
%     cprintf('red','WARNING: In Freq Domain: Ratio of ~%d --> less than 1 order of magnitude!\n\n',ceil(ErrFreq/ErrFreqRand));
% end 
% 
% %ploting the random signal and original signal
% subplot(3,3,7:9);
% plot(random_rec_vec,abs(random_rec_freq),orig_vec,abs(orig_freq));
% title('Random signal and Original signal - Freq domain');
% legend('Random Signal','Original Signal');
% xlabel('Frequency (Hz)');
% ylabel('Amplitude');
% xlim([-fmax,fmax]);

%% 
 
